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Abstract 

The  operational  implementation  of  the  ensemble  Kalman  filter  Composite  Clock  (CC 
developed  by  K.  R.  Brown),  is  the  subject  of  this  paper.  Although  the  mathematical  background 
of  the  CC  is  well-known,  its  autonomous  and  robust  operation  in  a  clock  laboratory  requires 
further  modifications  of  the  CC.  The  suitable  initialization  of  the  first  ensemble  estimate  is 
discussed.  Due  to  the  fact,  that  only  difference  measurements  are  available,  there  is  a 
dimensional  degree  of  freedom  in  the  first  as  in  the  following  ensemble  estimates. 
Furthermore,  to  guarantee  a  robust,  consistent,  and  autonomous  estimation  of  the  ensemble 
clocks,  a  consistency  check  is  described.  It  outputs  a  so-called  consistency  matrix  which 
describes  the  active  and  non-active  clocks.  Since  the  identified  active  clocks  can  change  from 
measurement  to  measurement,  the  corresponding  Kalman  filter  (KF)  parameters  are  adapted  in 
a  consistent  way.  Finally,  the  non-active  clocks  are  estimated  based  on  the  KF  output.  The 
performance  of  the  OCC  using  an  ensemble  of  five  DLR  laboratory  clocks,  measured  over  a 
period  of  1  year,  is  outlined.  The  clock  ensemble  consists  of  three  high-performance  cesium- 
clocks  (HP5071A),  an  active  H-maser  (CHI-75),  and  a  GPS-disciplined  Rb  clock.  Although 
miscellaneous  anomalies  like  phase  steps,  outliers,  and  clock  exclusions  arise,  the  OCC 
autonomously  processes  the  five  ensemble  clocks  and  establishes  a  robust  system  time.  Analysis 
results  of  this  ensemble  with  different  types  of  clocks  when  applying  the  OCC  are  presented. 


SYSTEM  TIME  CONCEPTS  FOR  GNSS 

Besides  orbits  and  other  system  parameters,  a  Global  Navigation  Satellite  System  (GNSS)  generates  and 
provides  satellite  clock  offsets  to  the  users.  The  clock  or  time  offsets  are  with  respect  to  (w.r.t.)  the 
GNSS  system  time  (ST).  It  is  an  integral  part  of  the  GNSS  and  is  a  key  element  of  the  GNSS 
performance. 

There  can  be  extracted  three  requirements  on  the  system  time: 

1 .  Stability  of  ST :  no  impact  on  satellite  clocks  for  any  x  value 

2.  Robustness  of  ST:  stability  is  guaranteed  at  any  time  t 

3.  Metrology  of  ST:  representation  of  UTC  and  long-term  performance. 

Regarding  the  ST  implementation  which  shall  provide  the  listed  requirements,  there  can  be  distinguished 
two  concepts. 
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The  first  ST  concept  is  based  on  a  master  clock.  The  system  time  is  established  by  a  highly  stable  atomic 
clock,  e.g.  an  active  hydrogen  maser.  Its  short-term  stability  of  about  2xl0'13@l  s  generally  fulfills  the 
first  requirement.  However,  to  meet  the  second  requirement,  “robustness,”  the  master  concept  has  to  deal 
with  master  clock  failures  like  outage,  time  or  frequency  steps,  or  outliers.  A  steered  backup  master  clock 
can  be  installed  to  solve  these  issues.  In  case  of  a  detected  master  clock  failure,  the  system  time  is 
established  by  the  backup  clock.  The  challenge  is  to  provide  such  additional  hardware  and  technology  to 
steer  the  backup  clock.  Similarly,  the  third  requirement  is  fulfilled  by  steering  the  master  clock  w.r.t. 
UTC. 

The  second  ST  concept  is  based  on  ensemble  time.  The  system  time  has  got  no  physical  representation 
and  is  computed  by  a  timescale  algorithm  ([1]),  which  computes  the  time  offsets  of  the  ensemble  clocks 
to  the  system  time.  The  system  time  is  basically  understandable  as  a  weighted  average  out  of  the 
ensemble  clocks.  The  stability  can  be  assumed  to  be  better  than  the  best  ensemble  clock.  Regarding  the 
second  requirement,  “robustness,”  the  system  time  accepts  loss  of  single  clocks,  which  provides  high 
reliability.  Similar  to  the  master  concept,  in  order  to  meet  the  metrology  requirement,  the  ensemble 
clocks  are  steered  to  UTC. 

ENSEMBLE  TIME:  COMPOSITE  CLOCK  (CC) 

The  paper  focuses  on  the  timescale  algorithm  composite  clock  ([2-4]).  It  is  a  Kalman  filter  (KF)  which 
models  each  ensemble  clock  by  three  states  (Appendices  A  and  B). 

The  Kalman  filter  is  a  recursive  method  which  estimates  the  ensemble  at  time  point  tk  using  the  previous 
ensemble  estimate  and  the  measurement  of  time  point  tk  with  t  \—tk—  tk_x . 

Assuming  the  ensemble  consists  out  of  N  clocks,  the  inputs  of  the  k-th  Kalman  filter  iteration  are 

•  Z(tk)  e  RN1  measurement  w.r.t.  measurement  reference  clock  (MRC) 

•  X  (tk  , )  e  R3N  previous  ensemble  estimate 

•  cq-  ^  G  r3Nx3n  previous  KF  covariance. 

Every  KF-iteration  can  be  separated  into  three  steps.  The  first  step  is  the  prediction  of  the 
ensemble  X~ (tk)  =  OX (tk_} )  .  Next,  the  measurement  residual  w.r.t.  MRC  is  calculated: 


Z(tk)  =  Z(tk)-HOX(tk_l) 

with  H  =  H#mrcP  (Appendix  B).  From  KF-theory  it  is  known  that  Z(tk)  is  a  White  and  Gaussian 
process  with  zero  mean  ([14]).  The  conditional  mean  of  the  ensemble  is  calculated  by 

X(tk)  =  OX(tk_l)  +  K(tk)Z(tk) 


The  corresponding  Kalman  gain  is 

K(tk)  =  C~(tk_x)HT(HC~(tk_x)HT  +R)~l 
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with  C  (tk_t)  =  OC(4_,)( Dr  +Q(t)  . 

Notice,  the  notation  X (tk )  =  (X#1  (tk), X#N  (tk))T  is  used.  X#,(tk)  e  R3 is  the  estimate  of  clock  #i. 
Besides  the  ensemble  estimate  X  (tk )  the  iteration  outputs  the  KF  covariance 


C(tk)  =  C-{tk_l)-K(tk)HC~(tk_l). 


In  detail,  K.  R.  Brown  developed  several  fundamental  properties  of  the  composite  clock  ([4]).  A 
fundamental  result  is  to  understand  the  ensemble  estimate  w.r.t.  to  an  implicitly  defined  system  time 


XST(tk)' 


R3 


X(tk)  =  X(tk)-HXST  {tk )  +  Nrep  (tk )  (1) 

with  H  :=(/3---/3)r  eR3Nx3.  N rep{tk)  &  R3N  is  called  representation  error  (to  the  implicit  system 
time).  The  covariance  of  the  representation  error  Crep  (tk )  :=  Cov(Nrep  (t, ))  is 

(tk )  =  C(tk )  -H{HtCx  (tk  )H)lHT 

([4]).  A  proof  given  by  Brown  shows  that  the  KF  estimates  X (tk )  are  unaffected  using  either  Crep(tk )  or 
('(/,  )  (transparent  variation).  The  advantage  of  Crep(tk )  compared  to  C(t, )  is  that  its  matrix  entries 
converge  which  corresponds  to  the  implicit  steady  state  of  the  KF  composite  clock.  However,  the  regular 
KF  covariance  C(tk) does  not  converge. 


OPERATIONAL  COMPOSITE  CLOCK  (OCC) 

Although  the  Kalman  filter  composite  clock  is  mathematically  specified,  its  application  in  a  real- 
world  system  results  in  several  additional  challenges  on  the  CC.  The  stability  of  atomic  clocks, 
which  is  modeled  by  the  q-values  (Appendix  A),  is  not  guaranteed  for  any  time  point.  Clocks 
can  be  corrupted  by  different  types  of  anomalies:  time,  frequency,  or  drift  steps.  Furthermore, 
the  measurements  can  be  disturbed  by  outliers.  None  of  these  anomalies  is  accurately  modeled 
by  the  CC  itself,  thus  corrupting  the  ensemble  estimate. 

The  operational  composite  clock  (OCC)  is  the  extended  version  of  the  CC,  which  deals  with 
operational  issues  of  its  real-world  application. 


OCC  MODULE:  INITIALIZATION 

At  the  moment  of  the  first  available  measurement  Z(tt) ,  no  previous  ensemble  estimate  X (t()) 
and  covariance  C(f0)  is  available. 
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Option  I:  No  Statistical  Information 
KF  theory  ([5-7])  recommends  setting 

X(t0):=E(X(t0))  and  C(t0):=Cov(X(t0)) . 


It  is  a  fundamental  property  of  clock  processes  that  the  mean  and  covariance  of  X(tk)  is  unknown  for 
any  time  point  tk.  An  option  is  to  assume  no  statistical  information  about  X(t0)  and  set  X(t0)  :=  0  and 
to  scale  the  corresponding  clock  covariance  with  a  factor  /  :=  (/:  l2  l3) : 


(  1111 

1  3  ,  /  _  1  5  /  „  1  _2  ,  /  „  1  4 


qM)  := 


20 


8 

hw  hq?> 


i  ^ 


\qxT  +  l2q2  t  +  l3q3  t  l2q2  z  +  l3q3  z  l3q3  z 


hq?>  2 T 
i3q3z 


Plugging  the  qt(z)  of  the  clocks  together,  the  according  ensemble  process  covariance  is  Qt(z)  and 
C(t0)  :=  Qj  (z) .  The  initialization  depends  on  the  choice  of  1  and  is  called  option  I. 


Option  II:  Put  First  Measurement  Residual  to  Zero 

Using  the  properties  of  the  measurement  residual  w.r.t.  MRC 

Z(t1)  =  Z(t1)-HQX(t0) 

a  further  initialization  method  can  be  derived.  In  an  optimal  designed  KF,  the  measurement  residual 
process  shall  be  zero  mean  ([14]);  thus,  a  reasonable  choice  is  to  put  X  (t0)  :=  ju  with 


0  =  Z(t1)-HOju.  (2) 

The  equation  (2)  is  affine  linear  and  the  rank  of  H ®  is 

rank(HO)  =  N-l 

(without  proof).  The  dimensional  freedom  of  the  solution  is 

dim(ker(H®))  =  3A^  -  rank(HQ>)  =  2N  +  l. 

An  ad-hoc  solution  of  equation  (2)  is  to  put  clock  #i  with  i  ^  MRC 

X #i  (t0 )  :=  :=  ( Z(.  (t0)  + a  0  ()f  and  X#MRC(t()):=  jU#MRC  :=(a  0  0)T 
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with  steer  parameter  “a.”  Since  the  frequency  and  drift  states  are  set  to  zero  and  the  measurement 
residual  is  zero,  there  holds 


xifd  =  ®X(t0)  +  K(tx)Z(t ,)  =  OX(t0)  =  X(t0) 
for  any  starting  covariance  C(t0)  . 


It  remains  a  suitable  choice  of  C(t0)  .  Here,  the  authors  are  not  sure.  The  applied  option  is  to  calculate 
offline  the  steady-state  covariance 


Crep(tJ-=\\mk^Crep(tk) 

by  successive  iteration  and  put 

C(t0):=  mCrep(tJ 


The  m  parameter  is  used  to  scale  the  initial  covariance. 

Option  III:  Put  First  and  Second  Measurement  Residuals  to  Zero 

In  option  II,  the  initialization  method  puts  the  frequency  states  to  zero.  The  method  could  run  into 
trouble  if  the  ensemble  includes  clocks  with  frequency  offset.  Typically,  RAFS  and  SPHM  have  a 
deterministic  drift,  which  leads  to  frequency  offsets  to  the  size  of  lxlO'12  or  even  more. 

In  order  to  properly  initialize  an  ensemble  including  clocks  with  deterministic  drift,  option  I  is  extended 
to  use  two  measurement  residuals.  Set  X  (/,,)  :=  u  with 


'z(ty 

'  HO  " 

l°J 

KZ(t2)j 

KHOOy 

(4) 


The  equation  (4)  is  again  affine  linear  and  the  rank  is 


rank 


HO 

kHOOj 


=  2(N  -1) 


(without  proof). 

As  a  result  the  dimensional  freedom  of  the  solution  is 


dim(ker 


HO 

kHOOj 


)  =  3N-2(N-l)  =  N-l  +  3. 
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The  authors  recommend  uniquely  solving  the  equation  (4)  by  fixing  the  drift  states  of  all  clocks  to  zero 
and  fixing  one  clock  without  frequency  offset,  e.g.  ju#MRC  :=  (0  0  0)r  ,  to  zero.  The  solution  is  called 

li  e  R  ' '  .  Again,  a  control  parameter  a  can  be  used  to  steer  the  time  offset  of  the  implicit  system  time. 
Setting  b#l  (a  0  0)T  for  every  i,  it  holds  that 


b  =  (b#1 


b#N)  eker 


f  //O  ^ 


and  jU  +b  solves  the  equation  (4).  The  first  and  second  ensemble  estimates  are 


X(t1)  =  0(Mu+b)  and  X(t2)  =  ®<S>(jUu+b). 

In  particular,  X1#MRC =  X1#MRC (t2)  =  a .  Assuming  the  fundamental  theorem  of  Brown,  the 

arguments  of  option  III  holds  and  the  implicit  system  time  is  steered  in  the  same  manner  as  presented  in 
the  example  of  option  II.  However,  the  steering  quality  depends  on  the  representation  error.  Analogous 
to  option  II,  the  initial  covariance  is  set  to 


C(t0):=  mCrep(tJ. 

Example:  Initializing  the  DLR  Clock  Ensemble 

The  initialization  options  are  compared  using  the  ensemble  of  five  atomic  clocks  operated  at  the  DLR 
timing  laboratory.  The  clock  ensemble  includes  two  high-performance  HP  cesium  clocks  (clock/cesium 
#1  and  #3),  one  standard  HP  cesium  clock  (clock/cesium  #2),  a  GPS -disciplined  rubidium  clock  (clock 
#4),  and  an  active  hydrogen  maser  (clock  #5/AHM)  of  KVARZ.  The  ensemble  clocks  are  measured  w.r.t. 
to  clock  #5  using  a  time -interval  counter. 

The  outlined  options  I,  II,  and  III  are  investigated  with  m  =  2  (option  II  and  III),  1  = 
(Ixl02°,lxl010,lxl010)  (option  I)  and  a  =  0  (option  I/II/III)  to  initialize  the  ensemble. 

Figure  2  and  3  show  the  2nd-state  estimates  of  the  cesium  #1  and  the  AHM  for  all  options.  In  Figure  1, 
the  2nd-state  estimates  are  biased  around  lxlO"13  for  option  I  and  III.  The  option  III  is  shifted  about 
0.5xl0'13. 

Similar,  the  2nd-state  estimates  of  the  AHM  are  biased  between  6xl0"14  and  10xl0'14  in  the  case  of  option 
I  and  III.  The  bias  of  option  II  is  around  lxlO'14. 

Assuming  that  the  clocks  are  free  of  a  frequency  offset,  option  II  consistently  estimates  the  clocks.  It  is 
the  most  promising  initialization  method  in  this  situation. 
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2nd  state  estimate  of  Cesium  #1 
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Figure  1 .  2nd-state  estimate  of  Cesium  #1 . 
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Figure  2.  2nd-state  estimate  of  AHM. 
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Mesasurement  residual  of  Cesium  #1  w.r.t.  AHM 
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Figure  3.  Measurement  residual  of  Cesium  #1  w.r.t.  AHM. 


Figure  3  shows  the  corresponding  measurement  residuals  of  Cesium  #1  w.r.t.  to  AHM.  By  visual  check, 
the  residuals  of  option  I  are  biased  till  the  day  two.  The  residuals  of  option  II  and  III  are  almost  identical 
and  unbiased  from  the  beginning. 


OCC  MODULE:  IDENTIFICATION  OF  NON-  AND  ACTIVE  CLOCKS 

The  k-th  iteration  processes  the  measurement  Z{tk)  to  compute  the  actual  ensemble  estimate  based  on  the 

former  estimate  X  (tk_1 ) .  Z(tk)  can  be  corrupted  by  outlier,  time  step,  non-available  measurement  .or 

frequency  steps.  None  of  these  anomalies  is  modelled  by  the  Kalman  filter;  thus,  a  regular  processing  of 
a  corrupted  measurement  likely  disturbs  the  ensemble  estimate.  It  is  the  fundamental  task  of  the 
consistency  check  to  detect  such  kinds  of  anomalies  in  order  to  provide  the  robustness  requirement  of  the 
OCC. 

Consistency  Check  of  Measurements 
The  statistics  of  the  measurement  residual  w.r.t.  MRC 


Z(tk)  =  Z(tk)-H®X(tk_l) 
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is  applied  to  do  a  consistency  check  of  the  residual  ([2]).  The  covariance  of  Z(tk  )  is 


HC  (4  )HT  +  R  =  HC;ep (4  )HT  +  R . 


A  hypothesis  test  of  each  component  of  Z(tk) is  performed  ([2]).  Clock  #i  is  called  active  w.r.t. 
MRC,  if 


Otherwise,  it  is  called  non-active  w.r.t.  MRC. 

It  remains  to  test  the  measurement  reference  clock.  The  authors  employ  an  ad-hoc  method.  In 
the  case  of  less  than  two  active  clocks  w.r.t.  MRC,  the  measurement  reference  clock  is  called 
non-active  and  active,  vice-versa.  By  definition,  a  clock  whose  measurement  is  not  available  is 
called  non-active,  too. 

Consistency  Check  of  Re-referenced  Measurements 

In  case  the  measurement  reference  clock  is  identified  as  non-active,  the  measurement  residual  is  not 
Z  (4 )  processed  and  an  alternate  reference  clock  is  determined. 

The  measurements  are  re-referenced  to  a  so-called  trial  reference  clock  #1: 

Z*  (4 )  :=  Z, (4 )  -  z,  (4 )  =  X?  (4  )  -  xrc (4 )  +  Vt (4 )  ■ -  (X”  (4 )  -  X[ #MRC (4 )  +  V,  (4 )) 


=  X1*(4)-X1«(4)+^(4)-Vl(4) 


•  Z”RC  Vk):=  -Z, (4 )  =  xrc (4 )  -  XT (4 )  ■ - V, (4 )  • 

By  adapting  the  measurement  matrix  to  Hm  :=  HtP  (Appendix  B),  a  hypothesis  test  of  the  measurement 
residual  w.r.t.  to  trial  reference  clock  #1 


Z#/  (4 )  :=  Zm  (4 )  -  H#l<t>X  (4_, ) 

is  performed: 

\z#l (tk )|M  <  A^H#,C:ep(tk)(Hm)T  +2R)u  . 

The  same  method  to  separate  non-active  and  active  clocks  is  executed  as  in  the  case  of  the  measurement 
reference  clock. 

Performing  a  consistency  check  w.r.t.  to  every  ensemble  clock,  a  so-called  consistency  matrix  is 
computed.  The  i-th  row  indicates  the  applied  trial  reference  clock  #i  and  the  j-th  column  indicates  the 
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measurement  residual  of  clock  #j.  The  matrix  entry  with  index  (i,j)  is  set  to  1,  if  the  consistency  check  of 
clock  #i  w.r.t.  to  trial  reference  #j  is  consistent,  and  zero  otherwise. 


OCC  MODULE:  KALMAN  FILTER  ADAPTATION 

The  consistency  matrix  is  used  to  determine  the  so  called  Kalman  filter  reference  clock  (KRC).  In  case 
the  entry  of  the  measurement  reference  clock  is  1,  it  is  used  as  KRC.  Otherwise,  a  clock  with  diagonal 
entry  “1”  is  selected  as  KRC.  The  corresponding  row  vector  determines  the  active  clocks  (having  “1” 
entries).  The  measurement  and  measurement  matrix  are  adapted  in  a  way  that  the  non-active  clocks  are 
excluded.  Furthermore,  the  measurement  matrix  is  adapted  w.r.t.  KRC.  In  case  no  valid  KRC  can  be 
identified,  the  ensemble  estimate  and  covariance  are  predicted. 

Since  the  number  of  active  clocks  M  can  be  less  than  N,  the  former  ensemble  estimate  X (tk_x)  and  KF 
covariance  Crep  (/),_, )  ,  which  are  of  dimension  3N  and  (3N,3N),  are  adopted.  The  non-active  clocks  are 
excluded  from  the  former  estimate  X (tk_x)  to  define  M(tk_x)  e  R3M  . 

The  covariance  adaptation  distinguishes  two  cases.  In  the  first  case,  the  number  of  active  clocks  at  tk 
compared  to  tk_i  decreases.  The  row  and  columns  of  the  non-active  clocks  are  excluded  from  the  former 

covariance  Crep(tk_x )  to  define  P(tk_x)  g  R3Mx3M  . 


In  the  second  case,  the  number  of  active  clocks  increases.  C  (tk_x)  is  reset  to  the  steady  of  all  clocks 
(offline  computation).  Afterwards,  the  rows  and  columns  of  the  non-active  clocks  are  excluded  to  define 
P(tk_ i)  g  r3Mx3M  The  authors  observe  by  simulation  that  if  the  number  of  active  clocks  increases,  the 

covariance  does  not  reach  the  steady  state  reusing  the  previous  C  (tk_x)  .  That  is  the  reason  to  reset 
Crep(tk_x)  to  the  steady  state  of  all  clocks  and  exclude  the  non-active  clocks.  The  reduced  covariance 
converges. 

The  actual  estimate  M(tk)  and  covariance P(tk ) ,  using  transparent  variation,  is  computed  based  on 
M(tk_ i)  and  P(tk_t) . 

It  remains  to  set  Crep(tk)  e  R3Nx3N  #  Consequently,  the  rows  and  columns  of  the  active  clocks  of  C  (tk ) 
are  set  to  the  corresponding  ones  of  P(tk ) .  The  non-active  clocks  are  set  to  the  corresponding  ones  of 

Crep  (4-i )  • 

Notice,  the  entries  Crep(tk )  of  the  non-active  clocks  are  modified,  if  the  number  of  active  clocks 
increases.  In  this  situation,  the  number  of  active  decreases  and  some  entries  are  unmodified. 


OCC  MODULE:  NON-ACTIVE  CLOCK  CALCULATION 

It  remains  a  procedure  to  calculate  the  non-active  clocks  w.r.t.  to  the  implicit  system  time.  The 
procedure  depends  on  the  reason  of  on  non-activeness. 
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NO  MEASUREMENT  AVAILABLE 

In  case  of  a  non-available  measurement,  the  non-active  clock  #i  is  predicted:  X#l  (tk  )  :=  </>X#l  (tk_x)  . 


Consistency  Failure 

In  case  the  clock  is  non-active  due  to  violating  the  consistency  check,  the  procedure  depends  on  the 
consistency  of  iteration  tk_i.  If  the  clock  is  previously  active,  the  non-active  clock  #i  is 

predicted:  X#l  (tk)  :=  (f>X#l  ( tk_x )  .  It  is  assumed  that  the  consistency  check  fails  because  of  an  outlier.  If  it 

was  previously  non-active,  the  time  offset  w.r.t.  the  implicit  system  time  is  calculated  using  XxKRC  (tk) 

and  Z*KRC(tk): 

z*KRC  0 tk )  +  x*KRC  0 tk )  =  xx#i  (tk )  -  xx#KRC  (tk  )+vi(tk)  +  (XX#KRC  (tk )  -  xf  (tk )  +  nx#krc  (tk )) 

=  x? (tk )  -  Xf r (tk )  +  Vi(tk)  +  N*krc (tk )  =:  Xx#i (tk ) . 

The  second  and  third  states  of  clock  #i  are  predicted: 

+  and  Xt\tk):=Xl\tk_x). 


OVERVIEW  OF  OPERATIONAL  COMPOSITE  CLOCK  (OCC) 

Figure  4  shows  the  work  flow  of  the  Operational  Composite  Clock. 


ROBUSTNESS  AND  STABILITY  OF  DLR’S  OCC 

The  DLR  clock  ensemble  of  5  atomic  clocks  (example:  initialization)  is  processed  using  the  OCC. 

Handling  of  Outliers  and  Time  Steps 

Figure  5  shows  the  measurement  of  Cesium  #3  w.r.t.  AHM.  It  is  corrupted  by  several  outliers,  and  a  time 
step  occurs  at  t  =  68.74d.  Focusing  on  the  measurement  residual  of  Cesium  #3  w.r.t.  AHM  (Figure  6),  the 
outliers  and  the  time  step  can  be  visually  identified  and  are  detected  by  the  OCC  consistency  check 
module. 
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Figure  4.  Work  flow  of  Operational  Composite  Clock. 
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Figure  5.  Measure¬ 
ment  of  Cesium  #3 
w.r.t.  AHM. 


x  10s  Measurement  residual  Cesium  3  w.r.t.  AHM 
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Figure  6.  Measurement  residual  of  Cesium  #3  w.r.t.  AHM. 
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In  case  of  the  outliers  the  estimate  of  the  Cesium  #3  corresponds  to  its  prediction  (Figure  7).  However,  to 
deal  with  the  time  step  at  t  =  68.74d,  several  steps  are  executed  by  the  OCC.  As  outlined,  the  consistency 
matrix  is  used  to  choose  a  valid  KF  reference  clock.  At  t  =  68.74d,  no  valid  KRC  can  be  identified;  thus, 
the  ensemble  is  predicted.  The  same  enters  the  following  two  iterations.  Analyzing  the  fourth  iteration 
after  t  =  68.74d,  the  diagonal  entry  of  the  AHM  is  still  zero  and  set  to  non-active.  The  remaining  diagonal 
entries  are  1  and  can  be  used  as  KF  reference  clock.  The  time  step  of  the  AHM  is  finally  present  in  every 
clock  measurement;  thus,  it  is  dropped  out  by  the  re-referencement  method.  It  remains  in  the  AHM 
measurement  residual  w.r.t.  to  the  trial  reference  clock  and,  thus,  the  AHM  is  set  to  non-active  w.r.t.  the 
trial  reference  clock. 

The  OCC  module  KF  parameter  adaptation  excludes  the  AHM  and  set  the  KF  reference.  The  remaining 
active  clocks  are  estimated.  Figure  4  shows  the  estimate  of  Cesium  #3  -  it  is  not  impacted  by  the  time 
step. 


Contrarily,  the  estimate  of  AHM  includes  the  time  step  of  the  AHM  (Figure  8). 
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Figure  8.  First  state  estimate  of  AHM. 
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ALLAN  DEVIATION  OF  TIME  LAB  CLOCKS 

Both  the  KF  estimates  as  well  as  the  measurements  are  corrupted  by  anomalies;  thus,  the  empirical  Allan 
deviation  of  the  data  is  impacted,  having  no  exclusion  method.  The  Dynamic  Allan  deviation  ([15])  can 
be  used  to  handle  this  issue  in  another  way.  At  each  time  point  tk,  the  Allan  deviation  is  calculated  using 
a  fixed  window  of  data.  Data  windows,  free  of  anomalies,  are  used  to  calculate  the  Allan  deviation  of  the 
clocks. 

The  time-interval  counter  is  specified  with  2*10_11r_1 .  Figure  9  compares  the  stabilities  of  the  Cesium 
#1  w.r.t.  AHM  or  the  implicit  system  time.  The  stability  of  the  Cesium  #1  estimate  and  the  measurement 
w.r.t.  AHM  are  almost  identical.  This  indicates  that  neither  the  AHM  nor  the  implicit  system  time  impact 

the  stability  of  Cesium  #1.  Its  performance  is  extrapolated  as  6*10-12r-0  5,  worse  than  its  hardware 
specification  of  5*10-12r-0  5 . 
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Figure  9.  Stabilities  of  high  performance  Cesium  #1 . 


Figure  10  shows  the  stability  of  the  AHM  estimate.  Its  stability  is  extrapolated  as  1*10  12  r  0  5 ,  worse 

than  its  specification  of  2*10-13r-0  5 .  This  is  reasoned  because  the  cesium  clocks  are  part  of  the  implicit 
system  time.  The  stability  of  the  weighted  cesium  clocks  is  worse  than  the  stability  of  the  AHM.  As  a 
result,  the  stability  of  the  AHM  estimate  is  characterized  by  the  weighted  cesium  clock  stability. 
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Figure  10.  Stability  of  AHM. 


CONCLUSION 

The  paper  describes  required  methods  to  operate  the  composite  clock  in  a  real-world  environment  (Time- 
lab,  GNSS).  The  statistic  properties  of  the  measurement  residual  process,  which  is  a  zero  mean  and 
Gaussian  process  (KF  theory),  are  utilized  to  develop  consistent  initialization  procedures.  The 
recommended  initialization  procedure  puts  the  ensemble  state  in  a  way  to  solve  the  measurement  residual 
equation  for  either  one  or  two  measurement  residuals.  Testing  the  procedure  with  a  DLR  clock  ensemble 
of  five  atomic  clocks,  the  most  promising  results  for  that  ensemble  are  achieved  by  fixing  the  frequency 
states  as  well  as  the  drift  states  to  zero.  However,  initializing  a  clock  ensemble  including  clocks  with 
frequency  offsets  it  is  recommended  to  solve  the  measurement  residual  equation  using  two  measurement 
residuals.  As  a  result,  frequency  estimates  of  the  clock  ensemble  are  provided. 

Based  on  the  measurement  residual  process,  a  hypothesis  test  is  described  which  splits  the  ensemble  in 
active  and  non-active  clocks.  The  hypothesis  test  works  with  any  reference  clock  as  long  as  it  is  part  of 
the  ensemble.  The  Kalman  filter  is  executed  with  the  active  clocks  to  calculate  estimates  of  them.  The 
non-active  clocks  are  calculated  using  an  additional  method.  The  robustness  of  the  OCC  is  validated 
using  measurements  out  of  the  DLR  clock  laboratory.  Outliers  along  with  time  steps  of  the  reference 
clock  are  detected  and  processed  by  the  OCC  without  disturbing  the  remaining  estimates. 
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APPENDICES 


A:  CLOCK  AND  ENSEMBLE  PROCESS 


The  process  of  clock  #i  is  defined  by  ([4,9,16,17]): 

Xm(tk)  =  0X#i(tk_l)  +  W#i(tk) 


with  (/)  = 


r\  r  0.5  r2A 
0  1  T 

0  0  1 


■  W  l(tk)  is  a  White  and  zero  mean  Gaussian  process  with 


Cov(W*'(t, ))  =  Q*'(t)  ■ 


f  1  3  1  5  1  2  1  4  1 

qxT  +  q2  T  +q3  T  q2  t  +q3  x 

3  2(J  2  5  o 

1  3  1  2 

*  <h  T  + 

*  *  q3T 


.  (qi,  q2,  q3)-value  is  the  so- 


V  2 

called  q-value.  Correspondingly,  the  ensemble  and  noise  processes  are  defined 

by X(tk)  =  (X#1(tk),....,X#N(tk)f  and  W(tk)  =  (W*\tk\....,W#N (tk))T .  W*\tk),....,W#N (tk) are 

stochastically  independent  and 


Q(T)  =  Cov(W(tk)). 


It  is  O  =  IN  ®<J). 


B:  MEASUREMENT  PROCESS 


Define  the  matrix 


and 


1 

0 

0  ••• 

•••  0> 

0 

0 

0  1 

0  0 

...  0 

rNx3N 

,0 

...  1 

0  ol 
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H  = 


r  e  -e  ^ 

Ll  Cr 


e.._,  -e.. 


-er 


R 


N-lxN 


VeN-e,-J 

with  ei  e  RN  i-th  unit  vector  and  measurement  reference  clock  #r.  The  measurement  process  is 
defined  by 


Z(tk)  =  HrPX(tk)+V(tk) 

with  V(tk)  =  (Vftk)  •••  Vv_l (tk )f  and  Cov(V (tk ))  =  R .  The  measurement  matrix  is  defined  by 
H  -  HrP . 
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